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Abstract 

We exactly reformulate the lattice CP(iV — 1) spin model on a D di- 
mensional torus as a loop model whose configurations correspond to the 
complete set of strong coupling graphs of the original system. A Monte 
Carlo algorithm is described and tested that samples the loop model with 
its configurations stored and manipulated as a linked list. Complete ab- 
sence of critical slowing down and correspondingly small errors are found at 
D = 2 for several observables including the mass gap. Using two different 
standard lattice actions universality is demonstrated in a finite size scaling 
study. The topological charge is identified in the loop model but not yet 
investigated numerically. 
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1 Introduction 



For a few simple lattice field theories Monte Carlo simulation algorithms are known 
that are virtually free of critical slowing down. In these cases we have an almost 
qualitatively improved control over approaching the universal continuum or scaling 
limit which is essential for particle as well as statistical physics applications. One 
important line of such developments started from the Swendsen Wang cluster al- 
gorithm [1] for Potts, and in particular Ising models, in arbitrary dimension. In [2] 
the present author has shown how via embedded Ising spins these techniques can 
be extended to O(N) invariant sigma models. This has led to numerous high preci- 
sion simulations close to criticality in the literature. Unfortunately efforts to widen 
the applicability to other theories have since ended in frustration in many cases. 
A well known example concerning the CP(iV — 1) model family is discussed in [3]. 
Some theoretical understanding of the restriction to O(iV) is given in [I]. Signifi- 
cant progress in the simulation of CP(N— 1) models has nevertheless been made in 
the sequel, for example with the non-recursive multigrid ('unigrid') method in [5], 
[6]. In [7], [8] an interesting but somewhat indirect method via quantum models 
and their reduction from higher dimension was presented. Nonetheless CP(iV — 1) 
systems seem to remain a challenging testing ground for hopefully more system- 
atically generalizable attempts to boost the efficiency of the numerical evaluation 
of lattice field theory. 

A completely different recent research programme is based on the proposal of 
a Monte Carlo summation of the strong coupling series (in a certain simple form) 
which replaces the sampling of lattice field configurations. For many models of 
interest this series converges for any given finite volume and choice of parameters, 
and the in general infinite set of strong coupling graphs may be considered as an 
equivalent non-perturbative representation of the original lattice model. 'Worm' 
algorithms generalizing those described in [9] were the essential tool to allow for 
an efficient simulation of this ensemble. In this formulation criticality means that 
large graphs are important, which are far too numerous for systematic evaluation. 
The question if a Monte Carlo procedure can efficiently sample a sufficient subset 
seems rather different from the problem of collectively updating long distance 
correlated fields. Little other means than numerical experiments are presently 
available to answer this distinct and open question. We here extend the recent 
series of papers [TO], [H], [12] to give an affirmative answer also for CP(iV — 1). 

In particular in [10] many details for the treatment of non-Abelian models 
have been developed that are similar here. Therefore we have to constantly refer 
to this work and the present paper could not really become self-contained without 
excessive repetition. On the other hand here further progress is made for the 
simulation method that can be used to also render the O(N) model simulations 
even more efficient that what was described in [ID] . 
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The present paper is organized as follows. In section 2. we formulate the 
CP(N — 1) model in several discretizations and derive the equivalent loop models. 
For their simulation an algorithm is described in section 3 and numerically applied 
in section 4. In section 5. the loop formulation is generalized to include the 9 
parameter coupled to the topological charge, but related numerical experiments 
are left to possible future publications. We end in section 6. with some brief 
conclusions. 

2 CP (AT — 1) model as a loop model 
2.1 Explicit gauge field formulation 

The CP(iV — 1) model is formulated with spins on lattice sites with values in a 
complex N — 1 dimensional projective space. They may be represented by one- 
dimensional projectors or by complex iV-component unit vectors <j>{x) G C N , \<p\ — 
1, whose phases are irrelevant. One of the standard lattice actions [T3] in use can 
be written as 



- S[(f), U} = (3 Yfofa ^)<P\x)<p(x + fi) + U~\x, /i)0 f (x + /t)0(x)]. (1) 



We here sum over all links of a hypercubic periodic lattice in D dimensions with 
extent L M in direction /x = 0, 1,. — 1. We use lattice units in which the L M 
are integer and V = Y\ is the number of lattice cells or sites. In addition to the 
spin field <fi a U(l) gauge field U(x,fi) is included which can absorb local phase 
transformations of (f)(x). There also is a global SU(iV) invariance with in the 
fundamental representation. In a first step we consider U as a given 'background' 
field over which we integrate later. We define the following partition function with 
two contracted adjoint composite insertions 



where d\i is the normalized invariant measure on the 2N — 1 dimensional sphere 
of 2N real component unit vectors made from the real and imaginary part of <fi. 
The generalized Pauli-Gell-Mann matrices A a with a = 1,2,..., iV 2 — 1 are a basis 
for traceless N x N matrices that obey the normalization condition 




z 



(2) 



tr(A a A b ) = 25' 



(3) 



Using the easily proven completeness relation 
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we may verify that our normalization is such that for coinciding u = v 



Y[u,u;U] = 2(1 - l/N) 



\\Mo(z)) 



-S[4>,U] 



(5) 



holds and thus the ordinary partition function emerges up to a known factor. 

The single site generating function immediately follows from the one of the 
O(N) model M 



<i/i(0)e- ; 



J2F[n;N)tfjT 



n=0 



with 



F[n; N] 



r(N) 



n\T(N + n) 

On each link xfi (y = x + £i) we double-expand 



(6) 



(7) 



J[Ucf,Hx)4>(y)+U- 1 cf,^y)4>(x)] 



k,k=0 



P 



k+k 



k\k\ 



U k - k (<f>\x)<f>(y)y(<f>l(y)cf>(x)y, ( 



k (rhU 



then promote the integers to link fields k(x, /x), k(x, /i) and also introduce the 
differences 

3n(x) = j(x, n) = k(x, n) - k(x, //). (9) 

Auxiliary fields 



d(z) = y~]k(z - fi, fi) + k(z, jj) + 5 Z)U + 5 ZjV 



(10) 



and 



(11) 



count the number of factors <ft and <fi at each site in an expansion term. The 
measure ([6]) is U(l) invariant such that nonzero contributions only result if d(z) = 
d(z) holds at all sites. This condition is equivalent to the vanishing divergence or 
flux conservation 

d(x)-d(x) = d;j,(x)=0 (12) 

where d* is the nearest neighbor backward derivative. 
The relevant integral over <f> now is 



Y[u,v;U] = J2U 



k,k X V L 



ok(x,fj.)+k(x,fi) 

k(x, fi)\k(x, //)! 



nw)) 



(13) 



Precisely as in [IQj this spin integral can now be performed by replacing all 0, 0^ by 
derivatives with respect to sources d/dj\ d/dj and then using (Q on all sites. All 
SU(iV) indices get contracted and the various terms correspond to strong coupling 
graphs on the lattice. On each link we draw k(x,fi) lines with arrows in the 
positive direction and fc(x,/z) lines with opposite arrows. At each site z we again 
imagine all possible ways to saturate all derivatives as a 'switchboard' that connects 
(contracts) all surrounding lines pairwise with each other, always an incoming to 
an outgoing one. This involves the weight F(d(z); N) from the measure at that 
site. The 'connectors' are regarded as 2-vertices joining pairs of lines at sites. 
In this way oriented closed loops arise around which the contractions in internal 
space contribute a factor N per loop to the total weight. Two (chains of) lines 
connect pairwise the four inserted spins at u and v, they end in four 1- vertices. 
These lines are open geometrically (unless u = v) and are saturated in internal 
space by the A matrices. Because these are traceless there is only one nonzero 
pairing of the four spins: a positively oriented line pointing from a w-spin 
to a f-spin \<j>{v)\ and the other one from v to u. The result of these contractions 
is a factor 2 (A 2 — 1). In analogy to the 0(A) model in [10] we call these lines 
active loops, distinguished as the wu-loop and the vu-loop, as opposed to the 
remaining passive loops. Either active loop is called trivial if for u — v it entirely 
lives on that site, i.e. it has no lines and 2-vertices. We now regard Y[u,v;U] 
as the sum over all possible loop graphs A 6 £ 2 of which each consists of many 
arbitrarily overlapping, intersecting and backtracking oriented closed loops plus 
the two active loops ending in 1-vertices at u and v. At this stage we consider 
k, k,u,v as functions of A. All graphs in £ 2 satisfy (fT2|) . After counting the 
multiplicity of the terms corresponding to each graph in the way discussed in [TU] 
we arrive at the remarkably simple form 

y[U) = C^p-\u-v)Y[u,v-U] = P~V-^MA]Ar |A| Y[u jM (14) 

U,V Xfl 



with the weight 
W[A] ' 



S[A] 



(x,p,)+k(x,fi) 

XjJ, 



^fWrh) Y ^-W-VV). (15) 



In the exponent |A| is the number of closed loops in the configuration A (including 
the two active ones in our convention). The factor S [A] is the symmetry factor of 
the graph introduced in the erratum to |10| . The strictly positive weight p with 
the normalization p(0) = 1 has been first discussed in |TT] and will be chosen to 
our convenience later. 
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In the standard CP(iV — 1) model one now integrates over all U(x, //) indepen- 
dently link by link with the normalized U(l) measure. This enforces the constraint 
j(x,fi) = in ffT4|) on all links. We call this subset of graphs C 2 C C 2 - The flux 
represented by the arrows now has to vanish identically on all links instead of just 
being conserved at the sites. Then the loop partition function becomes 

Z = I DUy[U] = P^( u ~ v)W[A]N w . (16) 
J Ae£ 2 

2.2 Quart ic action formulation 

A second popular action [13] for the CP(N — 1) model contains only the field 
with the action 

-Sj0] = 2/3 g ^|0t( x )0( x + / 2)|^ (17) 

Xfl 

If we start from S q , the same steps as above immediately lead to the locally flux- 
less (k = k) graphs £ 2 

Z q = ^2p- l (u-v)W q [A]N^ (18) 

A6£ 2 

with the weight, modified by the slightly different multiplicities, 

r(iv) 



W q [A] 



S[A) 



TT i^ll (19) 



For either action the relation between the two point correlation of the original 
spin model and the ensemble (1T61) or ([15]) is easy to establish. With double angle 
expectation values defined by [W — > W q for (|T8|) ] 



« ( A )» = \ E P'^ u ~ vWiA]N^O(A) (20) 

we find 



Z 



In particular the susceptibility 



X = \ E(0 t (O)A>(O)0t( x ) A >(x)) = {(tr[0(O)0t( O )0(x)0t( x )]) - 1 

X X 

(22) 

can be measured as 

N - 1 ({p(u - v)) 

x = N «Q> • (23) 
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2.3 Adjoint formulation 



The action S q may also be rewritten completely in terms of adjoint composite 
fields 

J a (x) = f (x)A a 0(x). (24) 

We introduce 

- S Q [0] = p a J c (x)J c (x + fx) (25) 

XfJ, 

but note that in the path integral we integrate over <f>(x) as before. There is the 
trivial relation, again a consequence of (j3J), 

Sa[4>]=SMk=^ + ^(3a. (26) 

Thus S q and S a just differ by a constant shift, irrelevant for any correlation. In 
spite of this the all-order strong coupling expansions are far from identical. An 
expansion in powers of (3 a using S a would involve unoriented lines as for the O(N) 
model. A major difference arises however from the measure that would be relevant 
in this case 

/oo 
dfi((f))e baja = A[m,n;N}(b a b a ) m (d abc b a b b b c ) n (27) 
m,n=0 

with a source b a . In the adjoint representation there is (for iV > 2) a second 
invariant totally symmetric tensor beside 5 a b, namely 

46c = \ tr(A a A b A c + A 6 A a A c ). (28) 

Thus this expansion is more complicated with both two point and three point 
vertices available to contribute. Only for iV = 2 d a b c vanishes and the adjoint 
formulation falls back to the standard lattice formulation of the 0(3 = iV 2 — 1) 
model as is well known. But even here S q yields a new expansion and the numerical 
reproduction of 0(3) results in the CP(1) formulation is non-trivial. 



2.4 Nienhuis Boltzmann factor 

In [TU] we have locally modified the original Boltzmann factor by truncating its 
expansion on each bond by k(x,fi) ^ k max . The original model is recovered for 
large k max but a particularly interesting case arises for fc max = 1. Such modifica- 
tions and the conjecture of universality still holding have been introduced into the 
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literature a long time ago [H], |15] . Truncating in the S q formulation amounts to 
the replacement 

e -sM ^ = II [ X + iWWfa + A)!'] • ( 29 ) 

XfJ, 

The adjoint translation reads 

M q \<P] = (1 + 2{3 q /N) DV + Pa J a (x) J a (x + £)] (30) 



with the adjoint Nienhuis coupling 



& = ~~~ • (31) 

For f3 q ^ we find ^ (3 a ^ AT/2. The criticality observed for 0(3) in [10] corre- 
sponds to large positive (3 a and N = 2. This is obviously not reached for positive 
f3 q but only as q /* {—N/2). We now turn to the possibilities of efficient Monte 
Carlo algorithms for the CP(iV — 1) loop model just outlined. We will however 
find this to be restricted to positive (3 q and hence cannot at present implement the 
Nienhuis formulation as for the 0(3) model. 

3 Simulation of the CP(N — 1) loop model 
3.1 Algorithm R for real N 

We have arrived at a representation of the CP(iV — 1) model in terms of loops 
representing arbitrary strong coupling graphs. These configurations can be param- 
eterized by linked lists as described in detail in [10] with only minor adjustments. 
One of the two possible orientations of each loop (for example column 2 of the list) 
is identified now with the physical orientation of lines of the present graphs. The 
flag in column 4 must now allow for three different values to distinguish between 
passive loops, the uv loop and the v u loop. 

We first develop an algorithm to simulate the CP(iV — 1) loop ensemble ( IT6|) . 
We define a number of separate update steps such that each of them fulfills detailed 
balance. They will then be iterated in some order as the final update procedure. 
The moves are all Metropolis proposals for which we quote a ratio q which yields the 
acceptance probability min(l,g). We encounter cases where our a priori proposal 
probabilities are not symmetric. To achieve detailed balance this needs to be 
compensated in q, which is then not just the ratio of the Boltzmann weights of the 
two configurations involved. This in particular applies to the inclusion of S [A] as 
discussed in the erratum to [TU] . 



S 



I. Extension and retraction step: With probability p ext = l/(l+r ext ) we propose 
an extension step, otherwise the retraction step. In the extension branch we 
next choose with equal probability one of the 2D directions to move u to the 
corresponding neighbor u with a concurrent extension of both active loops. 
The proposal is accepted according to the ratio 

= 2£>/3 2 r ext p(u - v) 

qext (N + d(u)){N + d{u) + l)p(u-vY 1 } 

In the retraction branch u is pulled back over one link along the active loops 
with the ratio 

(jy + d {u) - 1)(N + d(u) - 2) p(u - v) 
q " et ~ 2D^r cxt p(u-v) {66) 

if both active loops lead to the same neighbor u, otherwise no move is made. 

II. Re-route step: We here want to change the SU(iV) contraction or line con- 
nectivity structure at u. Such a move is only considered here if u ^ v and 
d{u) > 2 holds. 

i. We pick one of the 2-vertices at u and propose to swap between the line 
pointing out of this 2-vertex and the uv loop emerging from its initial 
1-vertex at u. For the q ratio we need to distinguish further sub-cases. 

a) The chosen 2-vertex belongs to a passive loop. The latter then effec- 
tively gets inserted into the uv active loop at u, |A| is reduced by one 
and we accept/reject with q=l/N. 

b) The chosen 2-vertex belongs to the uv loop itself, which self-intersects 
at u. Then a section is detached from it forming a new passive loop, 
|A| goes up by one, q = N. 

No move is made if the chosen 2-vertex belongs to the v u loop. 

ii. As i. but the role of of the two active loops switched. 

In figured] we try to graphically illustrate the elementary moves. We iterate these 
steps — and similar ones focussing on v instead of u — according to the scheme 

1 Iteration = (l u ll u l v ll v ) NxV/2 . (34) 

The ergodicity of this algorithm is shown as usual. One has to convince oneself 
that any graph can be built starting from the empty one by repeating the above 
moves. The empty graph is indeed the configuration from which we will start all 
simulations. 
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II 



XX 



Figure 1: Illustration of elementary update steps I (left) and II (right). Dotted 
lines indicate the continuation with other parts of a graph A. The crosses are at 
the insertion site u. 



For a simulation with the quartic action (fl7|) we just have to replace (132]) by 



, _ 4D(3 q r cxt (k(l) + 1) p{u-v) 
gext (N + d(u))(N + d{u) + l)p(u-v) 1 j 



and Q33J by 



, (N + d(u)-l)(N + d(u)-2)p(u-v) 
qret ~ 4Dp q r ext k(l) p(u-vY 1 ' 

In the last two formulae I denotes the link between u and u, i.e. I = (u, p) if 
u = u + p and I = (u, p) if u = u — p, with p, denoting a unit vector in the positive 
fi direction. 

In comparison with [10J a few changes can be noticed. The probability p cx t is 
new here. In [10J we have effectively chosen the special value p ext = 1/(2.D+1). We 
found the greater flexibility here useful to prevent acceptance rates from getting 
small. In all runs reported about below we have taken p ext = 0.3 and found 
unproblematic high acceptance rates. In addition there are no analogs of the steps 
Hi, Ilii and III of [10]. We found that they are not essential and also the O(N) 
algorithm may be simplified correspondingly. As a more technical change we have 
slightly modified the handling of the threefold linked list. As before two linkages 
allow to travel along loops in both directions while a third set of pointers allows 
to efficiently enumerate the vertices at a given lattice site. As we delete vertices 
(free list entries) in the retract step, we now immediately re-adjust all pointers 
including the site related ones. 

3.2 Algorithm I for integer TV 
3.2.1 Why an improved method? 

We have extensively run the algorithm just described and successfully reproduced 
numbers from the literature [TB] [5] [B] for small and intermediate lattices (L < 200 
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in D = 2, N = 4) with permille errors using only moderate PC time. As for the 
0(N) simulations in [10] we have found that integrated autocorrelation times in 
units of iterations typically stay below one in our timeseries of 10 6 iterations. 

There nevertheless is a problem that is exacerbated here compared to O(N). 
The moves described above seem to be superficially local with a fixed number of 
operations around u and v. For the reconnect II we have to know however if the 
2-vertex that has been randomly chosen is part of a passive or an active loop, 
which represents nonlocal information that is locally available to us in the flag 
entry of the list [10J. The price to pay for this is that whenever sections of loops 
are detached or eaten up, these have to re-flagged. It turns out that thus, close 
to the continuum limit on large lattices, most of the CPU time is spent travelling 
around loops resetting flags. This problem was already discussed in [ID] . It is 
more severe for CP(iV — 1) than for O(N) because the loops are more numerous 
and probably longer for comparable correlation lengths and lattice sizes. 

To diagnose the problem in more detail we have measured the loop-size distri- 
bution 

|A| 

(37) 



i=i 



where the expectation value 



« O(A )» = <^%» (38) 

refers to the closed 'vacuum' graphs only. Further we have introduced and imple- 
mented the decomposition of each such graph A into individual closed loops 

|A| 

A = U A * ( 39 ) 
i=i 

of lengths (number of link-lines) |A»|. In [17] . inspired by percolation and random 
walk theory, the asymptotic scaling form 

i n oc n- T e~ en (40) 

has been proposed introducing the loop tension 9 and a power correction exponent 
r. We found that such a fit describes reasonably well data that we have generated 
with S q at N = 4 for correlations lengths m _1 = 2 ... .18 and size mL « 10. Our 
aim here was not a high precision estimation of 9 and r. A clean assessment of 
their systematic errors would require quite some effort as it depends on the chosen 
fit window and their mutual correlation. We content ourselves at present with 
quoting that our analysis suggests that our data are roughly described by 

rwl.8, 9^m 2 /17. (41) 
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If we now pick a random 2-vertex we hit long loops more frequently and it will 
be part of a loop of length n with a probability proportional to n£ n . The cost to 
re-flag such a loop will hence on average be 



Because # _1 grows roughly proportional to the squared correlation length in the 
continuum limit ('fractal dimension two' ) simulations with the R algorithm slow 
down asymptotically in a way similar to ordinary local methods. It is to be noted 
that this section of the code (the re-flag while-loop) is very simple and thus domi- 
nates on larger lattices only, and that loop simulations still have great advantages. 
In any case we have found that the CPU time per site for R grew roughly by a fac- 
tor 60 as the correlation length and L are scaled up by a factor 7. Slower memory 
access on larger lattices may however also have entered here to some degree. 

3.2.2 Elimination of slowing down 

To proceed it was instructive to draw some analogies to Fortuin-Kastelyn cluster 
based algorithms for the g-state Potts model. In an early proposal Sweeny [18] has 
worked with bond variables only. The Boltzmann weight then contains a factor q Nc , 
where N c is the number of percolation clusters implied by the bond configuration. 
It resembles our weight iV' A ', for instance by allowing to analytically continue in 
q. During a local bond update nonlocal information is required: does the status 
of the single bond change N c or not? This leads to the same kind of slowing 
down. Sweeny has designed a system of hierarchical express pointers to accelerate 
the travel around loops (like express trains of the New York subway). One could 
consider such an improvement also for R here. 

A much simpler solution of the problem was however the one of Swendsen and 
Wang [lj. They keep spin and bond variables and update them in alternating 
order. One may actually view the cluster- wise assigned spins as just a device 
to stochastically dissolve the nonlocal weight into local steps. Transferring this 
technique to our problem at hand we insert, for example in ( JT5"j) . the representation 



Here the enlarged phase space now consists of graphs A where each loop Aj con- 
tained in it carries a label «j that is freely summed over N values. The update 
is extended now to such configurations including a choice for all «j. In the last 
(primed) sum we omit a assignments where the indices of the two active loops 




(42) 




(43) 
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coincide. The reason for this small extra twist will become clear soon. The loop 
partition function now reads for example 

Z=f DUy[U] = - v)W[A] (44) 

and expectation values in this ensemble are formed in the obvious way. 

The simulation of the ensemble including the a values requires only minimal 
changes. The flag entries of the linked list are rededicated to now store the flavor 
values a G {1, 2, . . . , N} of each vertex that is inherited from the loop to which it 
belongs. The extend/retract step I is completely unchanged, flavor is just passed 
on in extensions. The re-route step II (at fixed a) becomes even simpler now and 
proceeds as follows: 

IF If u 7^ v and d(u) > 2 holds we pick one of the 2-vertices at u. If it carries 
the same flavor as the initial 1-vertex of the uv loop we (always) swap in the 
way described before, otherwise nothing is done. Then the same is repeated 
for vu loop at u. 

Because of the restriction a uv ^ a vu it never happens that a re-route now leads to 
a line connecting for example <ft(u) with 4>(u), a contribution not included in our 
definition of the class £ 2 or C 2 designed such that u, v map out the adjoint corre- 
lation. The larger class would lead to a correlation with a singlet part decaying to 
a known constant instead of zero. It could be canceled but we would presumably 
be left with more and unnecessary noise. 

For ergodicity we also have to periodically update the flavor assignments after 
a certain number of the steps just discussed. We expect the cost for one iteration 
to remain O(V) if we do this only after O(V) extend/retract steps. One obvious 
option for these new steps would be to identify all loops A« in A (including uv and 
vu) and re-flavor them randomly. Some thought shows however that for ergodicity 
it is already sufficient to only randomly re-flavor the two active loops to one of the 
N(N — 1) pairs a uv ^ a vu . We call this type of re-flavor step now III. Some short 
experiments identified the following combination as quite efficient 

1 Iteration = [(IJIJJI„) y/2 III]*. (45) 

In particular our experiments have shown that it is not profitable to re-flavor all 
loops. Most important, it was quite pleasant to find that the errors of relevant 
observables grow only by 10% or so if they are accumulated during the same 
number of such I-Iterations replacing the much more costly R-Iterations. A more 
detailed discussion of autocorrelations and (the absence of) critical slowing down 
will follow in the next section. 
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4 Numerical tests 



We here report data on a number of simulations in D = 2 dimensions. The reason 
for this restriction is that physical interest in the CP(N — 1) models and hence 
available data seem to be concentrated on this dimensionality. We expect our 
new formulation and algorithms to be generalizable to other dimensions without 
problems. 

We have implemented the algorithms R and I as C codes. The graphs A are 
encoded into a linked list as described in detail in [10J. In C a very natural imple- 
mentation uses structures of pointers that are dynamically created and erased. We 
have found however that the more static storage handling described in [10] offers 
slight advantages in speed and have hence returned to it in the final version. 

4.1 Validation with the action S q in large volumes 

We first simulate the action S q for N = 4 to become able to compare (and find 
consistency) with data in [TB] [5] || on a lattice by lattice basis. These simula- 
tions are summarized in table [TJ After each extend/retract we have continuously 



p q 


L 


X 


6 


6 




K 


2.3 


20 


11.064(9) 


2.5736(12) 


2.6045(10) 


7.7 


2.59585(24) 


2.5 


32 


26.214(28) 


4.4491(27) 


4.5097(23) 


7.1 


3.07843(18) 


2.7 


64 


79.57(11) 


8.7813(76) 


8.9018(61) 


7.2 


3.57230(10) 


2.9 


128 


275.13(50) 


18.521(22) 


18.837(18) 


6.8 


4.03968(6) 


3.1 


256 


930.2(2.2) 


38.087(62) 


38.639(50) 


6.6 


4.48085(3) 


3.3 


512 


2998.6(8.2) 


74.80(16) 


75.54(13) 


6.8 


4.90816(2) 



Table 1: Results of simulations with the quartic action ( TT7|) that are immediately 
comparable to published data. 



recorded the contributions to (1231) and to the time slice correlations 

G(t) = ((p(u - v)[5 t „ + S tn })) (46) 

where the 5 functions are L-periodic and we enhance the statistics by summing 
over both directions for L = L\ = L. As discussed in [10] (see also [11]) we 
took p(x) proportional to the free lattice propagator with a mass M close to the 
one expected for the simulated lattice. From successive pairs of time slices we 
determine an effective mass by solving 

G(t + 1) cosh(m(t + 1 - L/2) , , ln . 

Wr^ ^ = \} ( , T l , ^m = m cS (t +12). 47 

G(t) cosh(m(t — L/2) 
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Following [16] correlation lengths = l/m e e(tk + 1/2) are determined self-consis- 
tently such that k£k £ [£fc,tfc + 1], i- e - roughly at separations £ and 2£. 
In the column 



we list the average number of lines per link in the simulated graphs. Actually there 
are K lines of either orientation. As in the 0(A) model (see [10]) K is a direct 
measure of the internal energy due to the identity 



Each row in table [T] derives from 10 7 iterations of the I algorithm. We have 
stored our data as 10 5 blocks from 100 successive iterations each. With these 
blocks an ordinary error analysis with the tool [19] was carried out. Between these 
measurements, on average separated by 100 iterations, hardly any autocorrelations 
are detectable. The errors of our errors are thus at a percent level and all digits 
given in the tables are significant. In addition, using multicore PCs, we always 
simulate between 8 and 32 independent replica and monitor for acceptable Q- values 



To extract the number of loops |A| in I-simulations one has to implement an 
additional observable to obtain this information. In our original R-simulations 
it is is however available 'for free' and was measured. We have found values 
|A|/V = 1.25 . . . .1.05 slowly falling as (3 rises, thus 0(1) loop per site. The 'his- 
toric' parameter sets in table [1] have led to physical sizes of about mL 7. We 
have found that this is too small to see a convincing mass-plateau at our present 
precision level. We have therefore repeated all runs on larger lattices mL 10, 
see table [2j In figure [2] the effective mass on the largest lattice now shows a satis- 
factory plateau. Beyond the steep initial decay, we show only every sixth effective 
mass value to not clutter up the plot. The size of the relative statistical errors of 
m e s is separation independent as expected The growth very close tot — L/2 
has kinematic reasons since at L/2 the correlation has a minimum for any mass. 
We then decided on these larger lattices to quote mass values from a fit to the 
correlation over a window t = L/4 to t = L/2. This is the horizontal line error- 
band in figure El We refer to [10] for a discussion of the fitting procedure which 
we took over unchanged. The lattices of table [2] represent a scaling set where the 
lattice spacing changes toward the continuum 'at fixed physics'. The weight p was 
formed with M = 10/ L. The optimal p is observable dependent. We repeated the 
run L = 380 with p = 1. The resulting error in x — 922.89(57) is about four times 
smaller while the error in £ = 38.089(33) went up by more than a factor two. The 
execution time per site of an I iteration between the first and the last line of table 
[2] still goes up by a (modest) factor 1.7, which is probably an effect of memory or 




(48) 



A = 2/3 g (|0 t (x)0(x + /i)| 2 >. 



(49) 



[19j. 
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Figure 2: Effective mass as a function of time slice separation at f3 q = 3.3 and 
L = 780. 



cache access with 8 replica (cores) sharing the memory. The total CPU time of 
the run L = 780 was about 330 hours on one double-quadcore Xeon (2.27GHz). 

We close on a warning side-remark. With all our observables being positive we 
originally thought that a single precision simulation (about 20 % faster) would be 
sufficient. We then found however that in plots like figure [2] the errorbars scattered 
around a smooth curve to a degree that seemed implausible. With the transition 
to double precision throughout (it was always used for the data analysis) this effect 
immediately went away. The mean value of the mass in single precision was still 
compatible with the new value while x was found to be different beyond errors. 

4.2 Autocorrelations 

In ordinary simulations the errors of observables are influenced by both the vari- 
ance of the quantity and by its integrated autocorrelation time. The latter depends 
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p 


L 


X 


f (L/4 L/2) 


K 


2.3 


26 


11.015(9) 


2.5889(7) 


2.59502(19) 


2.5 


44 


26.015(29) 


4.4561(13) 


3.07771(14) 


2.7 


88 


78.99(12) 


8.8256(28) 


3.57190(8) 


2.9 


188 


272.78(53) 


18.574(7) 


4.03950(4) 


3.1 


380 


926.4(2.4) 


38.068(15) 


4.48080(2) 


3.3 


780 


2960.2(8.6) 


74.619(32) 


4.90818(1) 



Table 2: Data from simulations of the CP (3) model with mL = 10. The mass 
m = l/£ is determined by a fit to the correlation. 



on the algorithm in use while the former is a property only of the ensemble to be 
simulated. This distinction tends to get blurred somewhat for our strong cou- 
pling simulations as already discussed in [TTJ. The point is that histograms like 
0(x) = ((6 X;U - V )) in the simplest case are measured continuously during the up- 
date. An extreme view for the R algorithm would be that after each microstep, 
for example in (|34|) . we in this way measure all 0{x) (mostly implicit zero 

contributions). Then the usual division holds and we could obtain r- mt in units of 
microsteps, but would need to handle lots of data for this purpose. In practice 
we block however 0(V) such measurements and then both the variance of the 
blocks and the residual autocorrelations between them depend on the underlying 
algorithm's 'decorrelation power'. 

With this explained we cite some autocorrelation times for the I algorithm. 
To that end we have repeated the simulations of table [2] with only 10 6 iterations 
but storing the contribution for each of them separately (in particular all G(t)). 
Then an ordinary analysis [19J yields autocorrelation times in units of iterations 
(computational complexity 0(VN), like 'sweeps'). Results are given in figure [3j 
The dotted lines just connect data points for the same observable. Here r intjm 
refers to the fitted masses. Values based on m e s(L/4) cannot be shown as they 
would fall on top of T int)X , i.e. be close to 1/2 for all lattices which means no 
autocorrelations in our definition of r int . We here see the complete absence of any 
growth of r int as the continuum limit is taken. The slowest modes couple to the 
total graph-size K which however even speeds up in the continuum limit. In any 
case, the most naive (and important) conclusion from table [2] is simply, that we 
see an approximately constant relative error in £ with costs only growing linearly 
with the number of sites. 
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Figure 3: Integrated autocorrelation times in units of I-iterations. 



4.3 Universality and finite size scaling 

We conducted another series of tests for the N = 3 model in a finite size scaling 
situation. To compare with data in [7] we here switch to the second moment 
definition of the mass m%. It is based on ratios of momentum space correlations 
which can be measured in the loop ensemble by 

G(p) oc ((p(u — v) cos(p ■ (u — f )))) (50) 

which follows from Fourier transforming (I2T1) . Table [3] contains our new data. 
We consider the step scaling function [20] z(2L) versus z(L) for 

z(L) = m 2 (L) x L (51) 

with pairs (L, 2L) at the same 0. For large enough L this graph is expected to 
reach a universal continuum curve. Rather then tuning the smaller lattices as 
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act. 




m 2 (L)L\ L=32 


m 2 (L)L\ L=m 


m 2 (L)L\ L=12S 




2.25 


2.3979(18) 


3.6582(22) 


7.2640(30) 


s q 


2.4 


1.9896(18) 


2.4133(19) 


3.6631(24) 


S q 


2.5 


1.8351(18) 


2.1080(19) 


2.6905(21) 


s q 


3.0 


1.4285(19) 


1.5286(20) 


1.6575(20) 


s 


4.0 


2.1138(17) 


2.6974(19) 


4.6179(26) 


s 


4.2 


1.9366(17) 


2.2812(19) 


3.1972(22) 


s 


4.5 


1.7499(18) 


1.9646(19) 


2.3433(20) 


s 


5.0 


1.5448(18) 


1.6747(19) 


1.8607(20) 



Table 3: Finite volume results for L = 32, 64, 128 with actions (IT7|) and ([T]). 



in [TD] we here give only a more qualitative 'curve-collapsing' demonstration. In 
figure H] we show the pairs contained in table [3] together with data0 from refs. [7] , 
[8] using their lattice pairs between (32, 64) and (104, 208) produced with standard 
simulations using S q at (3 q = 2.25, 2.5. We see a very convincing close to universal 
curve. Some remaining 'roughness' from lattice artefacts expected at a level L~ 2 
can just still be anticipated, for example around z(L) ~ 2.4. 

5 Topological charge in the loop model 

In two dimensions the gauge field U(x, fx) gives rise to an integer topological charge 
Q. It has a straight-forward definition on the lattice and we may hence extend the 
action by the 8 term by including the phase exp(i9Q) with the U integrations. 
We parameterize 

U(x,fi)=e lA ^ x \ A,(x)e(-n,n] (52) 
and define the field strength 

Ffw = d^A v - d u A^ e (-4tt, 4tt]. (53) 
An integer topological charge density is then extracted by setting 

F^ u = [Ffu,] + 27rn^, [F^] G (-vr, vr] (54) 
and the global charge is given by 

Q = J2n 01 (x). (55) 



1 I would like to thank the authors of [S] for sending their data and allowing me to reproduce 
them here. 
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Figure 4: Finite size scaling function of the CP (2) model with data from two 
different actions and spin as well as loop simulations. Errors (horizontal and 
vertical) are at most comparable to the symbol sizes. 



We next Fourier expand 

oo oo 

e i9[Foi] = H(6) h)e lh[Foi] = H (Q> h)e ihFm (56) 

h=—oo h=—oo 

with = 6/{2n) and 

= sinjg-/Qj = sirup . 
V ; (6-h)n 0/2 V ; 6-2nh K J 

We remark that H is precisely the kernel that appears in connection with the 
sampling theorem [2TJ. There a continuous time signal with compact support in 
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frequency space is reconstructed from its values at equidistant discrete times. In 
our case 9 is continuous and the critical Nyquist sampling frequency is one in our 
units. If 9 gets close to an integer k we have 



H(k + e;h) = 5, 



kh 



(5f 



in a distribution sense [first insert H into a sufficiently convergent sum, then take 
e — > 0]. If we now promote h to a field h(x), use Y2 X -^^ = an d rearrange the 
U(x, /i) from plaquettes into in a link-wise order we arrive at 



Y[H(9;h(x)) 



l[[u{x,n)]- 



-s,j, u d*h(x) 



(59) 



Xfl 



where e^ v is the antisymmetric tensor with e i — +!• 
The expression 



X^L h 



\[H(e-h{x)) 



n*. 



'j l i.{x),£ !lv d*h(x) 



(60) 



.<;/( 



is now well-prepared for its use in fll4p and leads to 

Z e =J2 P^( u - ^[A]iV |A| 0[A], 
AeZ 2 



(61) 



where G depends on the integer flux j(x,fi) which we here regard as a function 
of A. The set £2 is such that flux conservation ( Fl2|) holds, which is a necessary 
condition for the constraint in to have solutions. 

In [22] a very nice discussion is given for a closely related representation of 
cr-models which we adopt now. The field h(x) is really associated with plaquettes, 
or sites in the dual lattice. Its integer values in an h configuration may be viewed 
as the height of a tower above (or below) its plaquette, like a 'lego-plot' of the 
experimentalists. Then the euclidean plane is decomposed into domains of equal 
height. Due to the constraint in flBUl) the links along their boundaries carry nonva- 
nishing conserved flux j(x, fi) whose value and sign is determined by the difference 
between the domains heights on both sides. As noted in [22] this is a solid-on-solid 
model with additional interactions. In the limit 9 — > only h(x) = survives. At 
9 = it on the other hand there is a reflection symmetry h(x) — 1/2 -h- — (h(x) — 1/2). 

We leave it to a future investigation to decide if or at which (/3, 9) in spite of the 
sign oscillations in the weight H it is possible to numerically investigate topology 
in the loop version of the CP(A^ — 1) model. 
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6 Conclusions 



We have constructed a loop re-formulation of the CP(iV — 1) model which required 
just a rather straight-forward extension of the steps used in the O(N) model in [TU] . 
With only minor changes we could obtain the transcription for both the quartic 
lattice action and for the formulation involving an additional U(l) gauge field. For 
the numerical simulation an important refinement was necessary to avoid critical 
slowing down in the CP(iV — 1) case: the stochastic rather then exact treatment 
of the weight iV' A ' where |A| is the number of closed loops in an arbitrary strong 
coupling graph. The exact algorithm could still be used to simulate the theory 
also at non-integer N, if desired, although with accepting some critical slowing 
down. 

One of the main reasons to invent and study the CP (AT — 1) models was their 
analogy with Yang-Mills theory in four dimensions with regard to both asymptotic 
freedom and the existence of an integer winding number. As we have argued that 
the loop model is an exact representation of the original model the observables re- 
lated to topology should possess an image in the loop model. We have constructed 
the effect of the term in the action related to the 9 parameter. It leads however 
to the appearance of negative weights in the loop model. It is deferred to future 
work to find out in which parameter range, by including the signs in observables 
or by further re-formulation, numerical calculations are possible. 
Acknowledgments: I would like to thank Tomasz Korzec and Peter Weisz for 
discussions and Burkhard Bunk and Stefan Schaefer for advice with computing 
and C. Financial support of the DFG via SFB transregio 9 is acknowledged. 
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